function y = func_ns_nonmod(q, zeta, m, CONSTS)

    eps = CONSTS.eps;
    eta = CONSTS.eta;

    y_vec =  func_total_nonmod(q, zeta, m, CONSTS); 
%     y_vec(1) = -1i*Inf;
%     y_vec(2) = -1i*Inf;
    y_s_vec = func_1_nonmod_analytical(q, zeta, m, CONSTS) + ...
              sqrt(abs(eta)/eps)*func_2_nonmod_analytical((sqrt(abs(eta)/eps)*q), zeta, m, CONSTS);
    y = y_vec - (y_s_vec); 
    y(isnan(y)) = 0;
    
    if(false)
        figure; plot(q, imag(y_vec), 'b-', q, imag(y_s_vec), 'r-', q, imag(y), 'k'); title('Im(y)'); 
        legend('K', 'K^{(S)}', 'K^{(R)} = K - K^{(S)}');
        figure; plot(q, real(y_vec), 'b-', q, real(y_s_vec), 'r-', q, real(y), 'k'); title('Re(y)');
        legend('K', 'K^{(S)}', 'K^{(R)} = K - K^{(S)}');
    end

end

